# Problem 9
# =========
# A Pythagorean triplet is a set of three 
# natural numbers, a < b < c, for which,
#    a^2 + b^2 = c^2
#
# For example, 
#    3^2 + 4^2 = 9 + 16 = 25 = 5^2.
#
# There exists exactly one Pythagorean triplet
# for which a + b + c = 1000.
#
# Find the product abc.

# a * b % 12 == 0
# either a is odd, or b is odd
# a and b are coprimes
# c & 1 = 1

# LOOP thru successively increasing multiples of 12
#     calculate all factors
#     LOOP thru each pair of factors
#         IF a is odd and b is even
#             THEN IF a^2 + b^2 is a perfect square
#                 yield (a,b,c)
#             END-IF
#         END-IF
#     END-LOOP
#     calculate next multiple of 12
# END-LOOP
import math

"""
    Returns a list of pairs of factors for the given integer.
"""
def factors(num):
    return [(n,num/n) for n in range(2,int(math.ceil(math.sqrt(num)))) if num%n==0]

"""
    Returns the greatest common denominator for the given
    integers. If the return value is 1, that means that 
    a and b are coprime.
"""
def gcd(a, b):
    if a>b:
        return gcd(b, a)
    elif b%a==0:
        return a
    else:
        return gcd(b%a, a)

"""
    Returns the truthiness that the two given integers
    are the length of the legs of a Pythagorean triangle.
    Criteria is:
    a. a and b are coprime
    b. a and b are neither both odd nor both even
    c. (a * b) % 12 == 0
    d. a^2 + b^2 yields a number that has an integer sqrt
       that is odd
"""
def pythagorean_triple(a, b):
    from decimal import Decimal
    if gcd(a, b) != 1:
        raise "Not coprime"
    elif (a&1)==(b&1):
        raise "Not of differing parity"
    elif a*b % 12 != 0:
        raise "Product not a multiple of 12"
    else:
        s = math.sqrt(a**2 + b**2)
        if Decimal(str(s)).compare(Decimal(int(s))) == 0:
            return (a, b, int(s))
        else:
            raise "Not part of a pythagorean triplet"

# a ranges up to 332
# b ranges up to 333
# max(a)*max(b) will give us our upper bound!
MAX = 332*333/12

f1000 = []
for p in factors(1000):
    f1000.append(p[0])
    f1000.append(p[1])

pt = ()
for ab in range(2,MAX): 
    for p in factors(ab*12):
        try:
            trip = pythagorean_triple(p[0], p[1])
            s = sum(trip)
            if s in f1000:
                pt = trip
                break
        except:
            pass

m = 1000 / sum(pt) 
answer = map(lambda x:x*m, pt)
print "=> ", answer
print "-> ", answer[0]*answer[1]*answer[2]
